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ABSTRACT 

In  this  pup*  i d incut  t!.r>  • important  applications  of  a class  of  (parametric) 
linear  complement  >*'i t y problem  ir:sing  indep.  ndently  from  such  diverse  areas  as  port- 

folio  analysis,  si  tngiivering  and  graduation.  After  explaining  how  the 

complement  irity  pr  .*i  ■ ' . ■ r . t in  these  applications,  we  perform  some  analytical  com- 

parisons (i  as. h1  on  ; .it  i on  count-  and  storage  requirements)  of  several  existing  al- 
gorithms for  solving  this  cla- . of  conplc  i . ntarity  problems.  We  shall  also  present 
computational  result:  to  supjiort  the  analytical  comparisons.  Finally,  we  deduce  some 
conclusions.  about  the  general  performance  of  these  algorithms. 

SIGNIFICANCE  AND  EXPLANATION 

For  a given  i -vector  q and  n x n matrix  M,  the  linear  complementarity  problem, 
denoted  by  (q,M) , in  to  find  an  n-vcctor  x such  that 

q 4 Hx  > 0,  x^O  and  xT(q  + Mx)  = 0 . 

The  parametric  linear  complementarity  problem  consists  of  the  family  of  linear  complemen- 
tarity problems  <q  4 ).p,M)  where  p is  an  n-vector  and  X is  a scalar  parameter  ranging 
in  a certain  subinterval  of  the  real  line. 

Since  its  introduction  in  the  sixties,  the  linear  complementarity  problem  has  be- 
come an  extremely  important  subject  in  the  field  of  mathematical  programming.  It  has 
numerous  applications  in  many  different  areas. 
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In  this  paper , we  discuss  three  important  applications  of  a class  of  (parametric) 
linear  complementarity  problems  arising  independently  from  such  diver sr  areas  as 
portfolio  analysis,  structural  engineering  and  actuarial  graduation.  Our  purposes  are: 
(i)  to  exDlain  hov  complementarity  problems  emerge  in  those  appl icat ions,  (11)  to  com- 
pare, analytically  and  numerically,  the  computational  efficiencies  of  several  existing, 
"special -purpose"  and  "qeneral-purpose"  algorithm,  to  solve  thi  clar  . of  complement ai - 
ity  problems,  and  (iii)  to  draw  some  conclusions  on  the  suitability  of  these  algorithms 
applied  to  this  particular  class  of  problems. 

AMS(MOS)  Subject  Classifications  - 90C20,  90CS0,  65K05 

Key  Words  - Appl icat ions,  (parametric)  linear  complementarity,  quadratic  proqramminq , 
portfolio  analysis,  index  models,  structural  engineering,  graduation, 
difference-equation  method,  general-purpose  algorithms,  special-purpose 
algorithms,  comparison,  computational  results 
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ON  TOE  SOLUTION  Of  SOME  (PARAMETRIC)  LINEAR  COMPILE MENTARITY  EMBLEMS 
WITH  APPLICATIONS  TO  PORTFOLIO  ANALYSIS,  STRUCTURAL  ENGINEERING  AND  GRADUATION 

Jong-Shi  Pang*,  Ikuyo  Kaneko**  and  Wayne  P.  Hallman** 

1 . INTRODUCTION 

for  an  n-vector  q and  n by  n matrix  M,  we  denote  the  linear  complementarity 

problem  (LCP)  of  finding  an  n-vector  x satisfying  the  conditions 

T 

q i Nx  > 0,  x^O  and  x (q  + Mx)  =0 

by  the  pair  (q,M)  . The  parametric  linear  compl ementarity  problem  (PLCr) , denoted  by 
the  triple  (q,p,M),  consists  of  the  family  of  LCP's  (q  + Xp,M),  where  p is  an  n-vector 
and  X is  a scalar  parameter  ranging  in  a certain  subinterval  of  the  real  line. 

The  purposes  of  this  paper  are:  (i)  to  ex>  lain  some  important  applications  of  a 
class  of  LCP's  and  PLCP'a  chai  .i  t.-rized  Ly  a certain  form  of  the  matrix  M involved; 

(ii)  to  analyze  and  compare  computational  ef  f i cic  neicu  of  several  existing,  "special- 
purpose"  and  "genera  1-pur  p >.e"  algorithms  to  solve  thi;  class  of  complementarity  preb)  ems  ; 

(iii)  to  present  computational  results  and  (iv)  to  draw  some  conclusions  on  the  suita- 
bility of  these  algorithms  applied  to  this  particular  class  of  problems. 

In  recent  years,  important  applications  of  the  LCP  and  PLCP  have  emerged  in  two  en- 
tirely different  areas.  The  application  of  (parametric)  quadrat ic  programming  to  portfolio 
analysis  is  well-known.  Recently,  one  of  the  authors  (191  lias  developed  an  efficient  al- 
gorithm to  solve  a class  of  simplified  but  practical  portfolio  analysis  problems.  The  al- 
gorithm is  based  on  a PICT  formulation  of  the  problem  and  is  expected  to  be  considerably 
more  efficient  than  existing  solution  methods. 

The  other  application  is  in  the  area  of  structural  engineering.  Maier  (see  (12)  e.g.) 
has  formulated,  as  a LCP,  tht  problem  of  d termining  t tie  behavior  of  structures  belonging  to 
a certain  broad  class.  Since  then,  many  variations  and  extensions  of  this  LCP  approach 
have  been  developed.  In  particular,  one  of  the  authors  (10)  has  proposed  a reformulation 
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arid  algorithm  thereof  which  solves  the  problem  considerably  more  efficiently  than  the  original 
method  of  Maier's. 

| An  important  and  striking  fact  is  that  the  LCP ' s arising  in  these  two  applications 

are  defined  by  matrices  of  exactly  the  same  form: 

(1.1)  M » t i GGT 

where  £ is  an  n by  n diagonal  matrix  with  positive  diagonal  entries  and  G is  n by 
) m with  m less  than  n . This  implies , in  particular,  that  the  special-purpose  algorithm  de- 

veloped in  (19)  for  the  portfolio  analysis  application  can  be  used  to  solve  the  structural 
engineering  problems. 

In  this  paper,  we  shall  be  concerned  with  the  LCP  and  PLCP  which  involve  a matiix  of 
the  form  (1.1)  Although  this  cl  ;s  of  problems  is  certainly  quite  restricted,  there  seem 
to  be  more  applications  other  than  the  two  already  mentioned.  Kot  instance,  we  have  rt  - 
cently  ob.  rved  that  a certain  problem  in  the  area  of  curve  fitting  and  smoothing,  or 
graduation,  leads  to  a LCP  with  matrix  M of  the  form  (1.1). 

The  form  (1.1)  implies  that  the  matrix  M is  symmetric  and  positive  definite.  There- 
fore, the-  LCP  or  PLCP  with  such  M can  be  solved  by  most  well-known  LCP  algorithms: 

Notably,  Lomke's  algorithm  (11),  Col t le-Pant zig  (2)  and  Graves'  (6)  principal  pivoting 
algorithm  for  the  (unparameterized)  LCP  and  a parametric  version  of  Graves'  algorithm  (see 
(11)  for  the  parametric  LCP. 

By  taking  advantage  of  the  special  foim  (1.1)  of  M,  among  other  special  properties 
of  the  problem,  the  special-purpose  algorithms  developed  for  the  particular  class  of  comp- 
lementarity problems  are  expected  to  be  significantly  more  efficient  than  those  general- 
purpose  algorithms.  The  importance  of  investigating  such  special  methods  is  justified  by 
the  fact  this  particular  class  of  problems  has  applications  of  considerable  practical 
values  as  mentioned  a)K>ve. 
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2.  PORTFOLIO  ANALYSIS  APPLICATION 


2.1.  Problem  Formulation.  Starting  with  this  section,  we  present  the  promised  application 
of  the  class  of  L'  's  and  TLCP’s  with  matrices  of  the  form  (1.1).  The  first  of  such  ap- 
plications arisen  in  the  area  of  portfolio  analysis. 

The  standard  portfolio  analysis  problem  is  one  where  an  investor  allocater,  his  initial 
wealth,  taken  to  be  unity  with  no  loss  in  generality,  among  n risky  securities  available 
in  the  market.  Assuming  the  security  returns  are  stochastic,  Markowitz,  in  his  pioneering 
work  [161  described  a theory  postulating  that  rational  investors  should  select  a portfolio 
from  the  set  of  all  "feasible"  portfolios  which  offers  minimum  risk  (measured  by  variance) 
for  a given  level  of  expected  return,  and  maximum  expected  return  for  a given  level  of 
risk.  Such  a portfolio  is  raid  to  be  efficient . The  objective  of  the  portfolio  analysis 
problem  is  then  to  determine  the  set  of  efficient  portfolios. 

The  above  model  of  Markowitz  is  known  as  the  mean-variance  model.  Since  its  intro- 
duction in  the  fifties,  the  model  has  dominated  a great  deal  of  the  literature  in  portfo) io 
analysis.  As  the  model  is  so  well-known,  wc  assume  that  the  reader  has  some  basic  knowl- 
edge- about  it  and  therefore  shall  not  discuss  its  financial  and  economical  i mplicat ions , 
among  other  things.  Instead,  we  proceed  directly  to  the  mathematical  formulation  of  the 
model . 

T 

Let  a particular  portfolio  be  defined  by  a vector  x = (x  , . . . ,x  ) where  x.  rep- 
1 1 n l 

resents  the  proportion  of  the  investor's  wealth  invested  in  security  i . Let  the  security 

T 

return:,  In.  represented  by  a vector  of  random  variables  R = (R,  ,...,R  ) with  means 

1 n 

T 

U (u,...  ,ii  1 and  covariance  matrix  V - (o  .)  where  o..  = Cov(R.,R  ) . Then  the 

In  ij  l]  l g 

fiortfolio  analysis  problem  is  to  find  that  value  of  x that  will 


(2.1) 


minimize 
subject  to 

and 


1 T ,, 

- x Vx  - 


T 

C X = 


T 

0p  x 

1 0 

1 . 


Here,  e is  the  vector  of  one's,  a given  positive  vector  of  upper  bounds  on  the 

pt o;>ort  ions  to  be  invested  in  the  securities,  and  0 is  the  coefficient  of  risk  aversion. 
By  varying  0 from  0 to  ■»,  all  efficient  portfolios  can  bo  determined.  The  upper 


bounding  constraints  exist  because  of  legal,  personal  or  institutional  reasons. 


Soim  com- 


ponents  of  the  vector  a may  be  infinity. 

For  every  fixed  value  of  0,  problem  (2.1)  is  a (convex)  quadratic  program  which  can 
be  solved  by  a number  of  algorithms:  notably,  Markowitz'  critical  line  method  [17), 

Wolfe's  algorithm  [27],  bonne's  algorithm  [11)  and  Von  Hohenbalken' s algorittun  12b] . Both 
Markowitz'  and  Wolfe's  algoritluns  are  parametric  (i.e.  will  find  the  solution  for  every  0 ), 
whereas  Lemke's  and  Von  Hohenbalken's  are  not. 

Despite  the  fact  that  (2.1)  is  a rather  simple  (parametric)  quadratic  program,  few 
(if  any)  of  the  currently  available  quadratic  programming  algorithms  can  deal  with  problems 
with  a large  numbs  r of  sucuntie.  , an  important  reason  being  the  fact  that  typically  the 
covariance  matrix  '•  r highly  dense,  thus  causing  computer  storage  problems.  In  addition 
to  this  computational  limitation,  the  general  mean-variance  model  also  confronts  the  in- 
formational difficulty  >.f  data  collection  and  statistical  estimation  of  the  covariances  of 
the  security  returns  in  large-scale  applications. 

To  simplify  the  computational  and  informational  complexity  of  the  general  model, 

Sharpe  [ 2 -1 ) , [7M  proposed  some  simplified  models  for  portfolio  analysis.  These  are  the 

index  models  which  art  by  far  1 he  most  widely  used  formulation  of  the  portfolio  analysis 
problem.  In  a typical  m- index  model,  it  is  assumed  that  the  security  returns  are  deter- 
mined by  the  val';.-  f m market  indices.  Specifically, 

(2.2)  K = u.  + b I + ...  + b I + c , i = 1 n 

i l ill  imm  l 


where  I ‘s  are  the  market  indices  which  are  random  variables  with  mean  zero,  b.  ,'s 
3 x3 

are  constant  , and  i 's  are  i andom  variables  with  mean  zero  and  standard  deviation  c.'s  . 

i 1 

These  latter  random  variables  c.'s  are  used  as  measures  of  the  errors  of  approximating 

i 

the  actual  security  returns  lb's  by  taking  linear  combinations  of  the  market  indices. 

The  model  further  assume  • 


(A.  1 ) 

(A. 2) 


Cov  ( t—  ,t  ^ ) = 0 for  i / j 

Cov(t  . ,1  .)  = 0 for  all  i,j  . 


Under  these  two  assumptions , the  covariance  matrix  of  security  returns  is  given  by 

V = l + BCB1 


!• 

i 

t 

I 

i 

I 


(2.3) 
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I 

i 


2 

where  T is  the  n by  n diagonal  matrix  whose  i-th  diagonal  entry  is  o.. 


B *- 


bn  blm 


b , ...  b 

nl  rim 


is  the  n * m matrix  of  constants  b.  . h 

*3 


and  C is  the  m by  m covariance  matrix  of  the  market  indices.  It  would  be  convenient 
for  us  to  assume  that  the  oVs  are  positive.  Since  C is  symmetric  and  positive  semi- 


definite,  we  may  decompose  C into 


betting  G - BF,  we  may  write 


C = IT 


V = £ + GC 


which  is  of  the  form  (1.1). 

In  the  rest  of  this  subsection,  we  sketch  how  problem  (2.1)  can  be  solved  via  the 
PLCP  approach.  To  simplify  the  discussion,  we  assume  that  all  component!  of  the  vector  a 
of  upper  bounds  are  infinity  and  rtfi'i  the  reader  to  [lb]  for  the  general  case. 

The  Karush-Kuhn-Tucker  optimality  conditions  for  (2.1)  ar-  to  find  a vector  x and 
scalar  1 such  that 


(2.4a) 

u - 

-Xe  - 0p  + Vx  > 0, 

X > 0 

(2.4b) 

T 

U X 

= 0 

(2.4c) 

T 

C X 

= 1 . 

To  solve  (2.4),  one  would  normally  treat  X as  a variable.  However,  one  may  also 
treat  it  as  a parameter.  In  this  latter  situation,  (2.4a)  and  (2.4b)  become , for  evei y 
fixed  0,  the  defining  conditions  for  the  PLCP  (-0ti,-e,V)  with  X as  the  parameter. 

Now  if  we  solve  this  latter  PLCP  for  all  values  of  X and  if  we  can  determine  a suitable 

* — ♦ * 

value  X such  that  the  solution  x(X  ) to  the  LCF  (-X  o - 0p,V)  also  satisfies  (2.4c), 

* * 

then  x(X  ) and  X will  furnish  a solution  to  (2.4)  for  this  particular  value  of  0 
The  actual  application  of  the  PLCP  technique  to  solve  (2.4)  for  all  nonnegative 
values  of  0 is  divided  into  two  stages;  the  first  stage  being  the  determination  of  the 
initial  solution  corresponding  to  0=0  and  the  second  stage  being  the  solution  of  the 


problem  for  all  positive  0 . We  have  seen  how  the  first  stage  can  be  carried  out.  The 
second  stage  starts  after  the  initial  solution  is  obtained.  The  parameter  0 is  intro- 
duced back  in  the  problem,  with  its  value  to  be  increased.  The  problem  now  involves  two 
parameters  A and  0 . Nevertheless,  it  is  possible  to  eliminate  A (using  (2.4c))  and 
express  it  in  terms  of  0 and  the  values  of  the  basic  x-variables.  Therefore,  the  problem 
is  again  left  with  one  parameter  aid  the  solution  process  can  be  completed  by  treating  0 
as  the  only  parameter,  for  a more  detailed  discussion  of  the  method,  see  [19). 

The  P1.CP  approach  described  above  is  applicable  to  the  portfolio  analysis  problem 
(2.1)  with  a general,  nonsingular  covariance  matrix  V . In  this  paper,  we  shall  be  con- 
cerned only  with  the  index  models 

2.2  Compact  Inverse  Principal  Pivoting  Algorithm.  The  purpose  of  this  subsection  is  to 
describe  a special-purpose  algorithm  for  solving  the  class  of  LCP's  and  PLCP's  having 
matrices  of  the  form  (1.1)  The  algorithm  is  based  on  a modification  of  a parameterized 
version  of  Graves'  principal  pivoting  algorithm  (6)  for  solving  a PLCP  (q,p,M)  with  a 
(symmetric  and)  positive  definite  in.it  i ix  M . Since  the  parametric  Graves'  algorithm  is  also 
the  basis  for  the  "condensed  Graves'  algorithm"  to  be  described  in  the  next  section,  it  is 
useful  for  us  to  review  briefly  t h i . . basic  parametric  algoritlim. 

Consider  a PI.CP  (q,p,M)  with  a general  (symmetric  and)  positive  definite  K . Ac- 
cording to  a fundamental  theorem  in  the  theory  of  the  LCP  [22] , there  exists  a unique 
solution  xP)  to  (q  + Ap,M)  fo.  every  A . The  parametric  Graves'  algoritlim  will  de- 
termine the  solution  x(A)  for  all  A in  a finite  number  of  steps.  Suppose  for  conveni- 
ence, that,  a value  A is  obtained  such  that 

(2.5)  q + Ap  >_  0 for  all  A A 

This  implies,  in  particular,  that  x(A)  = 0 for  all  A A . Now  we  want  to  determine 
x(A)  for  A > A . To  achieve  this,  we  perform  a minimum  ratio  test  to  determine  how  much 
the  value  of  A can  be  increased  without  violating  condition  (2.5).  Assuming  nonde- 
generacy, we  either  obtain  a unique  index  where  the  minimum  ratio  is  attained,  or  arrive 
at  a nonnegative  p In  the  first  case,  we  perform  a principal  pivot  on  M,  obtaining  a 
principal  pivotal  transform  of  M and  new  vectors  q and  p . This  then  completes  an 
iteration  of  the  algorithm.  In  the  second  case,  a further  increase  of  A can  only 


increase  the  values  of  the  current  basic  variables;  therefore  the  algorithm  terminates. 

For  a more  detailed  description  of  the  algoritlim,  see  [1J. 

Two  remarks  are  important . Firstly,  if  one  is  interested  in  the  solution  x(A)  for 
A lying  in  some  bounded  interval  only,  one  should  incorporate  an  additional  stopping  rule 
so  that  the  algoritlim  will  also  terminate  when  the  right  end  of  the  interval  is  reached, 
in  the  case  where  this  latter  situation  occurs  before  the  arrival  of  a nonnegative  p . 
Secondly,  an  ordinary  LCP  (q,M)  with  a (symmetric  and)  positive  definite  M can  be  solved 
by  the  parametric  Graves’  algorithm  as  well.  To  do  this,  set  up  th<  PLCP  (f,q-f,M)  where 
f is  an  arbitrary  nonnegative  vector  such  that  f ^ q . Clearly  X - 0 satisfies  (2.5). 
When  X reaches  1,  the  problem  (q,M)  is  solved. 

Compared  to  the  compact  inverse  principal  pivoting  algorithm  to  be  presented  si.  . ly, 
the  parametric  Graver'  algorithm  is  more  general  in  the  sense  that  it  solves  any  (P)LCP 
as  long  as  the  matrix  rs  (symmetric  and)  positive  (semi-)  definite  (or  has  positive  pirinci- 
pal  min  ■ ■■  '.  Being  general  in  nature,  it  does  not  take  advantage  of  any  possible  special 
structure  of  the  matrix  M,  except  possibly  sparsity.  In  (19),  one  of  the  authors  lias 
proposed  a modified  version  of  the  algorithm  in  such  a way  that  any  special  structure  of 
M will  be  preserved  and  may  be  used  profitably  in  the  entire  solution  process.  Roughly 
speaking,  this  is  achieved  by  replacing  the  pivot  steps  in  the  original  version  by  up- 
dating the  index  set  of  basic  vatiablcs.  Here,  we  shall  not  describe  the  modified  algorithm. 
Instead,  we  present  as  Algorithm  1 below  its  specialization  to  a PLCP  (q,p,M)  with  M 
of  the  form  (1.1).  We  need  to  introduce  some  notations  first.  Let  A be  an  n by  m 
matrix  and  a {l,...,n},  p c (l,...,m}.  By  A „ we  denote  that  submatrix  of  A con- 
sisting  of  entries  whose  row  and  column  indices  belong  to  a and  6 respectively.  If 

a (l,...,n)  (S  = {l,...,m}),  we  write  A „ (A  ) to  mean  A . Similarly,  if  q 

. P a.  ap 

is  an  n-vcctor,  by  q^  we  mean  that  subvector  of  q whose  components  have  indices  be- 
longing to  a . We  attach  subscripts  to  the  identity  matrix  I to  denote  its  dimension. 

We  assume  condition  (2.5)  holds  for  some  X . 

Algorithm  I.  Compact  inverse  principal  pivoting  algorithm. 

Step  0.  bet  X-  min  (-q./p.)  and  let  k be  a minimizinq  index.  Let  a = (l,...,n) 

‘ old  l i 

Pi  < 0 

\{k)  and  P = (k)  . Go  to  Step  1. 
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5teP  1 • Solve  the  systems  of  linear  equations  for  o & and  o 

(2.6) 


(Im+  ‘V^BB  <V>  <°  t’' 


qB  pe 


and  compute* 

(2.7) 


((vTewv  ig«.,TiJv 


qB  rrt 

s = q - Go  , t = p - Go 


— - P If  *a  - 0 and  i °>  set  * (*>  = 0 and  xg(A)  = -l~g  (sg  + Xt  ) for  all 

* — \,yd  and  terminate.  Otherwise  determine  the  new  critical  value 

Ei  s ■ 

<2'8)  Xnew  Iain{min  i f a,  t.  < 0},  min{-  ~ : i c 6,  t > 0))  , 

i i 1 

and  let  tin  new  critical  index  k be  a minimizing  index.  Put  x (A)  = 0 and  x (A)  = 

-1  OB 

^BB^B  + XV  t0r  X 6 1Aold'  Xnew>  ’ Set  Xold  = Xnc„'  qo  to  Step  3‘ 

ste£_3.  Replace  a and  B by  a u{k)  and  B\(k}  (a\{k)  and  Bu{k})  respectively, 

depending  on  k c 6(a)  . 


The  stopping  criterion  in  Step  2 is  set  up  for  solving  the  PLCP  completely.  Other 
stopping  rules  should  of  course  be  included  to  deal  with  other  situations. 

The  larges t poi  t ion  of  work  required  in  Algoritiun  1 is  the  execution  of  step  1.  In 
fact,  the  overall  efficiency  of  the  algoritiun  is  completely  determined  by  how  fast  Step  1 
can  be  earned  out.  in  what  follows,  we  derm  strate  how  this  important  step  can  be  executed 
most  efficiently.  To  this  end,  let 


(2.9) 


A(B)  = Xm+  <Gg  / ^ (G6_) 


and  write 


(2.10) 


s(B> 


t (B)  p 


The  matrix  A ( 6)  is  symmetric,  positive  definite  and  l.as  spectral  ladius  at  least  1 . 
Moreover,  A(B)  is  of  order  m which  is  a constant. 
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r 


(2.12) 


r8.k  _ B,k 

f,  = On 


8,k 


with  n tht  (unique)  solution  to 
(2.13) 


P V T 

A(B)  n (G  ) 

k . 


(:-.(£<)  ), 


(2.14) 


(t(P)), 


K(P,k)  = 


. : , k 

1 + f 


K,(e'k)  • -T—BTk 

<l°k  'k 


with 


(2 . IS) 


-) 


if  k e B ' \ B 
if  k c B \ 8 1 ; 


and  finally,  o is  the  k th  diagonal  entry  of  E . 


Formula  (2.31)  shows  that  if  s(B)  and  t(B)  are  known,  one  can  compute  s(B') 


and  t(B')  by  first  solving  for  n^’k  in  (2.13),  and  then  computing  f,  'k,  K(£,k)  and 


B,k 


K'(B,k)  using  (2.12),  (2.14)  and  (2.15).  A fast  approach  to  solve  for  n ' is  to  apply 
Cholesky  factorization  to  the  matrix  A (8)  . Observing  that 


(2.16) 


A (6  • ) = A(B)  + (Do”2(Gk  )T  (Gk  ) 


which  identifies  A(B’)  as  a rank-one  modification  of  A(B) , we  may  apply  a very  efficient 
algorithm  developed  in  (5)  to  update  the  Cholesky  factors  of  A(B’)  We  rephrase  this 
algorithm  in  Subroutine  I below. 


Assume  that  the  Cholesky  factors  of  A(B)  are  known,  viz  A(B)  = LDL  where  L 


is  a unit  lower  triangular  matrix  and  D is  a positive  diagonal  matrix,  we  wish  to  de- 

---t 


termine  the  factors  A(B')  = LDI,  In  the  subroutine  below,  the  matrices  L and  D 
have  overwritten  L and  D respectively. 
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Subroutine  X. 


(Gill  et  al . 15))  Updating  the  Cholesky  factors. 


1. 

LX* fine  gamma 

= toe,  abd 

k 

2. 

For  j = 1 , . . 

. ,m,  compute 

dsave  = d(j) 

temp  = gamma  * .’(j) 

d ( j ) = d(j)  4 temp  * Z(j) 

beta  temp/d(j) 

gamma  dsave  * gamma/d (j) 

Z(k)  = Z(k)  - Z(j)  * i(k,j) 


k = j + 1 , . . . , m 


i(k,j)  = f(k,j)  4 beta  * Z(k)  J 

Since  both  A(B)  and  A(B')  have  spectral  radii  at  least  1,  the  above  procedure 

is  numerically  stable  Once  the  Cholesky  factors  of  A(B)  are  known,  the  solution  for 

B,k 

H in  (2.13)  involves  the  solution  of  two  unit  triangular  systems  of  linear  equations 

B k 

and  n divisions.  More  precisely,  the  computation  of  n ' consists  of  the  following 
three  steps: 


Subroutine  II  Solving  (2.13). 

T 

Step  1.  Solve  for  5 in  the  system  1*6  = (G^  ) 

Step  2.  Compute  D ^5  . 

B k T B k 

Step  3.  Solve  for  n * in  the  system  L n ' = 6 

We  now  describe  how  the  above  analysis  is  incorporated  in  tin.  implementation  of 

Algorithm  1.  Initially,  8 is  the  singleton  (k);  s(8)  arid  t(B)  are  computed  from 

4 k T 

(2.11)  by  setting  B'  = 8 and  8 = 4>  . Note  that  s($)  = q,t(4)  - p and  n = (G^  ) 

In  general,  at  the  beginning  of  each  iteration,  we  have  s(B)  and  t(B)  available.  The 

algorithm  then  goes  through  the  ratio  test  to  determine  if  termination  is  reached  or  a 

new  critical  index  k is  found.  In  the  latter  case,  a new  index  set  B‘  is  obtained  by 

inserting  or  deleting  the  index  k to  or  from  B . We  then  use  Subroutine  I to  update  the 

B k 

Cholesky  factors  of  A(8),  use  Subroutine  II  to  solve  for  n ' and  finally,  compute  the 
new  s(B’)  and  t(B')  by  (2.11).  This  completes  an  iteration  of  the  algorithm. 
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2.3.  Analytical  Comparison  of  Algorithms.  In  this  subsection,  we  make  some  comparisons 
between  the  compact  inverse  and  some  general-purpose  I,CP  algorithms  to  solve  the  (parametric) 
LCP  having  a matrix  of  the  form  (1.1).  As  "representatives"  of  general-purpose  algorithms, 
we  have  chosen  Lemke's  algoritlim  for  the  (nonparametr ic)  LCP  and  the  parametric  Graves' 
algor ithm  for  the  PLCP.  Clearly,  the  compact  inverse  algoritlim  is  applicable  to  both  the 
LCP  and  PLCP. 

As  criteria  of  the  comparison,  we  shall  consider  the  computer  storage  requirement 
(for  the  input  data  and  "basis  inverse  representation")  and  the  number  of  operations  (multi- 
plications and  divisions)  required  to  solve  the  problem.  Considering  the  fact  that  all 
three  algorithms  are  based  on  essentially  the  same  principles,  i.e.  simple*  pivoting,  re- 
sults with  respect  to  these  two  criteria  should  provide  a reasonable  measure  for  the  over- 
all performance  of  the  algorithms. 

Needless  to  say,  the  performance  of  an  algorithm  critically  depends  on  how  it  is 
implem.  nt  •••! . For  the  purpose  of  the  comparison,  no  consideration  is  given  on  sparsity  of 
problem  data.  Indeed,  in  portfolio  analysis  and  the  engineering  application  (53), 
the  data  are  highly  dense.  Also,  for  the  sake  of  simplicity,  we  assume  that  no  basis  re- 
inversion will  be  performed  in  the  course  of  calculations. 

A point  in  which  major  differences  could  occur  is  the  representation  of  basis  in- 
verse. In  the  compact  inverse  algoritlim,  it  is  assumed  to  be  exactly  as  specified  above. 

In  Lemke's  and  Graves'  algorithms,  we  assume  the  standard  product  form  using  an  "eta-file", 
as  described  in  14,  p.  264). 

The  comparison  results  are  summarized  in  Table  2.1. 


r 


compact  inverse 

IiOmke 

Gi aves 

; — 

storage 

3 

input  data 

q.p.l.G:  n(m+3) 

* 2 

q+X  p , M : n +n 

q,p,M:  n^  + 2n 

basis  inverse  , 
4 

represent  at i on 

L:  m(m-l)/2 

D:  m 

„ . 5 

eta- file  : T'n 

6 

eta- file  Tn 

t ot  al 

n(mt  3)+m(m+l) /2 

2 

n + n (T  + 1 ) 

2 

n + n(T+2) 

operation  counts 

. th  . 

k — iterat  ion 

(after  k-1 
pivots) 

Subroutine  I : m“+4ni 

2 

Subroutine  II:  m 

s and  t : n (m+2) +2 

ratio  test;  n 

ratio- test:  n 

eta-vector:  n 

HUS : n 

pivot  column:  kn 

ratio-test:  n 

pivot  column:  (k-l)n 

eta-vector : n 

RHS:  2n 

tot  cil 

O 

T [n (m+3)  t-2  (m+1)  ] 

h”  er-if)  , , 

nr^  _r 11 

n 

T(T-  1 ) 

2 

+3Tn 

Notes: 


Table  2.1  Storage  and  operation  counts  for  three  algorithms 

* 

1.  Lemkt  's  algoritlun  is  applied  to  the  problem  with  A = A 

2.  Only  real  arrays  are  taken  into  account. 

3.  Symmetry  of  M could  be  used  to  reduce  the  storage  requirement  in  l*emkc's 
and  Graves'  algorithms. 

4.  In  practice,  the  length  of  the  eta-file  can  be  controlled  by  reinverrion 
and/or  other  devices. 

5.  T'  - total  number  of  pivots  in  I«emke's  algor itlim. 

6.  T - total  number  of  pivot  r.  in  G aves'  a Igor  itlim. 

7.  The  operation  counts  for  the  ratio-tests  are  upper  bounds. 

From  the  al>ove  comparison  results,  the  following  two  observations  can  be  made: 

(i)  The  efficiencies  of  the  two  general -purpose  algorithms  are  about  the  same 
provided  that  the  numbers  of  pivots  to  solve  the  problem  for  each  of  the  two  algorithms 
are  close  to  each  othei * (there  is  no  reason  to  believe  that  one  algorithm  tenr-  ites  with 
consistently  fewer  or  more  pivots  than  the  other). 


1 The  number  of  pivots  in  Lemke's  algorithm  depends  on  the  choice  of  the  "coverinn  vector”. 
If  the  covering  vector  is  chosen  to  be  q ♦ Xp  which  is  nonnogativo  then  T = T'  holds. 
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(ii)  The  relative  performance  of  the  special-  and  general-purpose  algorithms  depend 
on  the  values  of  T (or  T')  and  m/n;  the  advantage  of  the  special-purpose  algorithm  be- 
comes more  significant  as  m/n  becomes  smaller  and  T larger. 

2.4.  Computational  Results.  In  this  subsection,  we  report  the  computational  results  with 
the  solution  of  some  variations  of  the  portfolio  analysis  problem  (2.1)  under  the  assump- 
tions of  an  m-index  model.  Before  ['resenting  the  results,  we  point  out  some  further 

T 

characteristics  of  the  matrix  ),  + GG  appearing  in  this  particular  application.  Fust  of 
all,  the  matrix  G tends  to  be  highly  dense  because  there  is  no  reason  to  assume  that  the 
indices  are  uncorrelated.  Secondly,  the  ratio  m/n  is  typically  very  small.  In  a poten- 
tial large-scale  application,  one  might  expect  n = 1000  and  m = 10  . Thirdly  and  lastly 
the  nuittber  n could  be  very  large  in  practice. 

The  experiments  performed  in  the  present  and  latci  t ion  wore  all  done  on  the 
UN1VAC  1110  at  the  Madison  Academic  Computing  Center  at  the  University  of  Wrsconsin- 
Madison.  All  the  input  data  were  generated  r a ; '.only.  The  c input  at  io  . in  this  section 
were  done  in  double  precision. 

The  main  purpose  of  the  experiments  below  was  to  test  the  capacity  and  efficiency  of 
the  compact  inverse  algorithm  in  solving  practical,  large-scale  portfolio  analysis  problems 
The  fact  that  the  problem  data  tend  to  be  highly  dense  makes  the  application  of  general- 
purpose  algorithms  impractical,  if  not  impossible,  for  solving  large-scale  problems  since 
in  these  algorithms,  all  nonzero  elements  of  M must  be  stored.  For  this  and  other 
reasons,  we  shall  not  report  here  any  computational  results  using  the  general -purpose  al- 
gorithms. See  conclusion  (ii)  in  Section  2.5  and  footnote  2 below. 

Three  sets  of  experiments  were  performed,  each  of  which  corresponded  to  solving  a 
variation  of  the  portfolio  analysis  problem.  A number  of  problems  corresponding  to  dif- 
ferent values  of  n were  solved  in  each  set  of  experiments.  For  each  value  of  n,  the 
integer  m was  allowed  to  vary  from  five  to  thirty,  in  the  scale  of  five,  in  order  to 
test  its  influence  on  the  efficiency  of  the  special-purpose  algorithm. 

The  first  set  of  experiments  was  designed  to  solve  the  following  problem: 

IT  T T 

(2.17)  minimize  ^ x (!)  + GG  )x  + q x 

subject  to  x > 0 

T 

and  e x = 1 


which  corresponds  to  problem  (2.1)  for  a fixed  value  of  9,  with  V of  the  form  (1.1) 
and  no  explicit  upper  bounds  (i.e.  a " for  all  i).  The  compact  inverse  alqorithm 
was  used  to  solve  the  PLCP  formulation  of  (2.17).  Two  representative  sets  of  results, 
corresponding  to  two  values  of  n,  are  contained  in  Figure  2.1.  We  have  used  time/iter- 
ation  as  a measure  of  the  efficiency  of  the  method.  The  timing  is  exclusive  of  input/out- 
put. 

The  second  set  of  experiments  was  concerned  with  solving  the  pr  Mem  (2.17)  in  para- 
metric form,  i.e. 

1 T T T T 

(2  -18)  minimize  — x (X+GGJx+qxt  Op  x 

subject  to  x > 0 

T 

and  e x - 1 


The  compact  inverse  algorithm  was  used  to  solve  tin  PLCP  formulation  of  the  problem.  One 
representative  set  of  results  is  presented  in  Figur,  2.1.  In  computing  the  time/iteration 
in  these  experiements,  only  the  • . ond  stage  of  the  PLCP  approach  is  taken  into  considera- 

tion . 

The  third  set  of  experiments  war,  designed  to  solve  problem  (2.18)  with  explicit  upper 
bounds : 


1 T _ 

T.  T 

minimize 

2 x <-  + 

GG  ) x 4 q x 4 

subject  to 

a > 

x _>  0 

and 

T 

e x 

= 1 . 

The  portfolio  analysis  problem  (2.1)  has  q -■  0 . Nevertheless,  q was  not  set  to  be 
zero  in  the  experiments.  An  extension  of  the  compact  inverse  algorithm  (19]  was  used  in 
this  instance.  The  upper  bounds  were  all  set  equal  to  1/20  . Two  representative 

sets  of  results  are  shown  in  Figure  2.2.  The  timing  is  computed  as  in  the  last  set  of 
experiments. 

Several  important  points  about  the  above  three  sets  of  experiments  are  summarized 

below: 

(i)  In  all  these  experiments,  the  values  of  m were  set  not  to  exceed  thirty  in 
order  to  be  consistent  with  the  smallness  of  the  ratio  m/n  in  most  practical  situations. 

(ii)  In  solving  (2.18)  (and  (2.19)  as  well)  by  the  PLCP  approach,  the  tiroe/iteration 


in  the  lir-t  tat  slightly  less  than  the  time/itoj  at  ion  in  the  second  stage.  This  is 
b.  aute  tti  tht  ond  stage,  somewhat  luor . quantities  have  to  be  computed. 

(lii)  In  addition  to  those  problems  whose  results  are  plotted  in  Figures  2.1  and 
/.2,  we  hive  solved  ■ ver.il  other  problems  of  the  same  nature.  All  of  them  produce  curves 
simii.n  t<>  those  m the  two  fiyures.  Therefore  they  are  omitted. 

(iv)  Tn-  r<  .nits  of  the  experiments  are  consi  .tent  with  the  operation  counts  of 
the  compact  inverse  algorithm  given  in  Table  2.1.  In  particular,  the  table  indicates  that 
the  operation  count  per  iteration  of  the  algorithm  is  independent  of  the  total  number  of 
iterations  and  quadr.it  io  in  m . These  two  facts  are  substantiated  by  the  two  figures 
(and  other  results  lot  reported  hero). 

( v)  The  total  numbers  of  iterations  required  to  solve  the  problems  (2.17)— (2.1°) 
depend  very  much  on  the  exact  type  and  sice  of  the  problems.  The  largest  number  of  itera- 
tion' (edit . Iis  1 in  the  experiment:.)  for  (2.17)  is  72  when  n - 300,  m = 30;  that  for 
(2. IB)  i-  ill  (excluding  tht  first  stage)  when  n BOO , m = 30;  and  that  for  (2.39)  is 
201  (excluding  the  first  stage)  when  n BOO,  m = 20  . 

2.B.  Conclusions.  Based  on  the  analy:  i.  and  computational  experience  we  have  gathered, 
we  may  draw  the  following  conclusions  about  the  suitability  of  algorithms  to  solve  practica' 
large-scale  portfolio  analysis  problems  of  the  kind'  given  by  (2 . 17) - (2 . 19)  . 

(i)  The  comp. let  inverse  algorithm  is  capable  of  handling  problems  of  potentially 
very  large  size  with  a minimum  storage  space  requirement  and  in  a reasonably  fast  manner. 

(ii)  Under  the  assumptions  that  n is  large  and  m/n  is  small,  which  is  the  case 
in  many  practical  situations,  it  is  evident  that  the  compact  inverse  algorithm  is  much 
preferred  to  the  general-purpose  algorithms.^ 


2 We  did  perform  some  numeric.,]  experiments  comparing  computation  times  required  by  the 
compact  inverse  and  I,emkp*s  algorittims  on  a number  of  problems  similar  to  (2.17).  We 
chose  not  to  present  the  results  here  since  the  superiority  of  the  former  is  clear  from 
the  operation  counts  analyst ' in  Table  2.1,  and  since  the  reliability  of  such  analysis 
had  been  established  by  the  computational  results  shown  earlier. 
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STRUCTURAL  KNUIKELKING  APPUC.MWK 


As  mentioned  in  Section  1,  the  (parametric)  LCP  with  a matrix  of  the  form  (1.1)  ap- 
pears in  the  nonlinear  structural  analysis  initiated  by  Maier  (see  e.g.  {12]).  Maier's 
L/ZP  model  war.  later  1 efonnulated  by  one  of  the  authors  in  terms  of  what  he  calls  a n by  dn 
LCP. 

In  this  section,  we  first  outline  how  a fundamental  problem  in  structural  engineering 
is  formulated  a;,  the  complementarity  problems.  After  describing  the  principal  pivoting 
algorithm  develoj  d in  {8]  to  deal  with  the  n by  dn  LCT,  we  jhall  present  some  results 
on  comparisons  of  algorithms  to  solve  the  structural  problem;  including  those  of  a system- 
atic computat  ional  expei  intent . 

3.1.  The  SI  r u ' urciJL  Problem  and  Its  Formulation.  The  formal  at  ion  outlined  below  is 
essentially  due  to  Maier.  To  minimize  technicality,  we  shall  treat  a simplified  model: 
we  refer  to  Maier  112],  and  Kaneko  {9,10]  for  the  details  and  for  more  general  approaches. 

A fundamental  problem  in  structural  engineering  is  to  determine  the  behavior  of  a 
structure  subjected  t.o  a set  of  loads.  It  is  assumed  that  the  behavior  of  the  entire- 
structure  can  be  represented  by  the  values  of  stresses  and  strains  dt  certain  (finitely 
many)  points  in  th*  structure.^  We  shall  call  them  critical  points  and  let  n be  the 
number  of  criti  al  points  present  in  the  structure  under  consideration. 

For  simplicity,  it  is  further  assumed  that  the  stress  and  strain  at  each  of  the 
critical  point-,  are  one-dimensicmal . 

The  methods  we  shall  deal  with  are  designed  for  a certain  particular,  but  practically 
impor  tant  < l aj  tructur*  railed  (piecewise  11  neai  » e j as  * plasti  In  partic  ul&r , 

reinforced  concrete  frames  belong  to  this  class. 

Under  appropriate  conditions,  the  set  of  mechanical  principles  which  governs  the 
behavior  of  the  clast  i c-plastic  structure  can  be  represented  by  the  following  (nonlinear) 
complement ai  i ty  problem- 

(3.1a)  q - NsE  t h(x)  - NZNx  » 0 

3 Those  points  may  represent  cross  sections  in  a l>oam,  or  finite  elements  in  a continuum, 
and  so  on. 
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(3.1b) 


x > 0 


Y E 

(3.1c)  x (q  - Ns  -*  h(x)  - NZNx)  = 0 

Here,  q is  a nonnecjative  n-vector,  N is  a n by  n diagonal  matrix  w*.;.  plus  or 

E E 

minus  ones  as  diagonal  elements,  s is  an  n-voctor  such  that  Ns  is  positive,  h is  a 

T 

function  of  the  form  h(x)  (h t (x, ),..., h (x  ))  where  h.(x  ) is  a concave , strictly 

11  n n 3 3 

monotone  increasing  real-valued  function  of  x ^ with  h_.(0)  - 0,  j = l,...,n,  and 

T 

finally  Z is  an  n by  n matrix  of  the  form  V.  -BSB  , where  S and  B are  m by 
m and  n by  m,  respectively. 

In  (3.1),  all  data  can  he  determined  by  mechanical  properties';  of  the  structure  under 
consideration  and  the  loads  aj plied  to  it.  In  particular,  m is  th2  number  of  critical 

points  with  a certain  property  (called  r-  dund.tr. * ) and  S denotes  the  stiffness  matrix 

with  respect  to  th<  redund«»nci<  thus  n < n and  S is  positive  definite. 

Suppo:  x i:.  a solution  of  (3.1).  From  it,  the  vectors  of  strtsr.es,  s,  and 

strains,  v,  are  obtained  by 

(3.2)  s = sE  + ZNx 
and 

(3.3)  v = Cs  + Nx  , 

where  C is  the  n by  n f K xibil i t y matrix  (with  respect  to  all  critical  points) . 

In  the  following,  wo  shall  attempt  to  explain  — very  briefly  — some  mechanical 
principle?  which  give  rise  to  (3.1). 

Letting  q’  - q + h(x)  and  using  the  expression  (3.2),  we  may  consider  (3.1a)  as 
requiring  the  "stress  point",  s,  lie  in  th»  convex  polyhedron 

(3.4)  E --  {s:  q’  - Ns  >_  0}  . 

I 

I 

The  set  E,  called  the  elast i • doma 1 n , defines  the  reqion  such  that  the  "plastic  activity" 
takes  place  only  if  the  stress  i>oint  reaches  the  boundary  of  E . j 

The  plastic  activity  ir.  represented  by  x,  the  vector  of  plastic  mu)  t ipliors , and 
affects  the  structural  behavior  in  two  ways.  First,  it  contributes  to  the  values  of 
stresses  and  strains  directly  through  (3.2)  and  (3.3),  respectively.  The  value  of  x 
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also  affects  the  value  of  q,  hence  the  shape  of  the  elastic  domain,  through  h . Ttu 
map  h is  referred  to  as  a (work*)  hardening  mle. 

The  plastic  flow  rules  specify  conditions  x must  satisfy;  ttasr  ai. 

(3.5)  x >_  0 

(3.6)  x1(q'  - Ns)  = 0 . 

The  boundary  hyperplanes  of  E represent  yielding  limits.  Thus,  the  rules  (3.5)  and  (3.(>) 
(togethei  with  s e E)  state  that  the  j11'  plastic  multiplier,  x , can  take  on  a josi- 
tive  value  only  if  the  stress  point  reaches  the  jth  yield  limit  (and  yielding  occurs 
there).  The  conditions  (3.5)- (3.6)  and  that  s c E give  (3.1) 

Finally,  a remark  on  a parametric  problem.  In  the  representation  (3.1),  the  effect 
of  the  applied  loads  is  accounted  for  in  the  form  of  the  vector  sK,  which  is  called  that 
°f  linear-elast i cresponses  to  the  loads.  It  is  often  of  importance  to  determine  the  com- 
plete evolution  of  stresses  and  strains  dining  a loading  process  where  the  vector  of  loads 
of  the  foim  If  is  applied  at  time  \ for  \ e [0,  A]  with  \ > 0 . This  can  be  done 
by  solving  the  parametric  form  of  (3.1),  given  as  follows,  for  all  X in  (0,  X)  : 

(3.7a)  q - XNsh  1 h(x)  - NZNx  >_  0 

(3.7b)  x > 0 

(3.7e)  x1  (q  - An.sK  + h (x)  - NZNx)  - 0 . 


In  this  section,  we  shall  only  considei , for  simplicity,  the  nonparametric  problem  (3.1), 
except  foi  the  description  of  the  solution  procedure. 


3.2.  LCP  Formulations.  By  assuming  that  h.(x^)  is  piecewise  lineal,  j = n , 

we  can  linearize  the  problem  (3.]),  in  two  alternative  ways,  so  that  pivoting  methods 
can  be  applied  to  solve  it. 

lirst  of  all  we  note  that  if  h(x)  = Hx  for  some  (positively  diagonal)  matrix  H , 
i.e.  if  every  fc  is  linear  then  (3.1)  itself  is  a LCP  with  the  matrix  11  - NBSBTN,  which 
can  be  put  in  the  form  (1.1)  (a  factorization  of  S must  be  obtained) 

In  the  following,  we  assume  that  tv(x^)  has  d (>  1)  linear  pieces  as  depicted  in 
Figure  3.0  for  j = l,...,n  . Such  h^(x^)  is  specified  by  giving  values  of  u*  > 0 , 

1 = 1 , . . . ,d  and  r > 0,  i = l,...,d-l 


L 
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Under  the  above 
that  x solve:  (3.1) 
Maier : 


assumption  on  h, 
if  and  only  if  a 


it  can  be  verified  (by  elementary  manipulations) 
dn-vector  x solves  the  following  LCP  due  to 


(3.8a) 


T E - T T 
q - P Ns  + (H  - P NBSB  NP)x  0 


(3.8b) 


X > 0 


(3.8c) 


xT(q  - PTNsE  + (H  “ PTNBSBTNP)x)  = 0 , 


with  the  relationship  x = 

,-f  dnxdn  . , 

H c R , are  specified 


Px 

as 


. In  this  problem, 
fol lows : 


the  new  data, 


q c R 


dn 


P € K 


n*dn 


and 


qd(j-l)+l  q3 

k k 

qd(i-l)+k+l  qd(j-l)+k  + Uj  rj' 


j = 1 n 


k = 1 d-1 


P = 


H. . =0, 
*3 


i j 


d(j-l)+l,  d(j-l)+l 


= u . >0 


d(j-l) +ktl,d ( j- 1 ) +k+l 


k k+1 
u . u . 

-r-3 .—  > 0 , k = 1 , . , 

k k-1 

u . - u . 

3 3 


j = ’ 1 , n 


. ,d-l 


Again,  after  factorizing  S we  can  put  the  matrix  in  (3.8)  into  the  form  (1.1).  Wc 

would  like  to  emphasize  here  that  by  passing  from  (3.1)  to  (3.8),  the  size  of  the  problem 

has  been  increased  by  a factor  of  d . 

The  linearization  scheme  employed  by  one  of  the  authors  [10]  is  based  on  a represen- 
tation of  the  piecewise  linear  function  h given  by 

d . . 

(3-9)  h (x ) = [ U1  x1  , 

i=l 
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with  i.  ,.di  t 1 'll 


(3.10) 

(3.11) 

and 


, i i.T  iH 
(r  - x ) x = 0, 


i = 1, . . . ,d-l,  x >0 

i = 1 , . . . ,d-l 


dn  i.cr 


x = l X1  , 
i*=l 


when  x t a , r (i 


..i1)1,  and  U1  diag{u* , . . . ,u* } , i = l,...,d  . 
n In 


Sul  at  1 1 , it  in  (3.1)  by  (3.9)  along  with  the  side  conditions  yields  the  n by 

4 


(3.12a) 


(3.12b) 


q - NsK  + 1 (U1  - NZN)x3  > 0 

i=l 

0 < x’<  r’,  i = l,...,d-l,  xd  > 0 


(3.12c) 


(3.1 2d) 


tJ  ^ 

xJ (q  - NsK  3 l (U1  - NZNlx1)  = 0 
i=l 

i i T i+1 

(r  - x ) x =0,  i = 1,..  ,,d-l 


Since  -Z  is  j os i ' ive  semi  -definite  and  0 is  diagonal  with  positive  diagonal 

elements  , i 1,  . .,d,  it  follows  that  the  n by  dn  matrix 

1 d 

(U  - NZN,  . . . , U - NZN) 

has  the  .-j^-propi  i ty  (See  Kaneko  (8))  and  thus  (3.12)  has  a unique  solution,  which  can 
be  computed  by  the  algoritlim  developed  in  (8)  (described  in  the  next  subsection). 

The  obviotr  advantage  of  the  n by  dn  formulation  (3.12)  is  that  it  involves  only 
n linear  constraint  , except  for  the  simple  upper  bounds  on  some  of  the  variables  which 
can  be  treated  as  in  the  standard  upper  bounding  technique  in  linear  programming. 

3.3  The  Conden-'-d  Gniv:,'  Algorithm  In  the  following,  we  shall  describe  the  principal 
pivoting  algorithm  developed  by  Kaneko  (8)  to  solve  an  n by  dn  LCP  with  a matrix 


4 Note  that  under  the  conditions  (3  10)  (3.11),  the  complementarity  (3.1c)  is  equivalent 
to  (3.  12c) . 

Id,  . 

f>  Tlie  fact  that  u.  > ...  > u.  > 0,  ] - l,...,n,  follows  from  the  assumed  monotonicity 

and  concavity  of1  h (x  ) 

3 3 
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having  the  ^-property.  The  algorithm  was  called  the  condensed  Craves*  algorithm  for  the 
reason  explained  in  [8J . The  algorithm  given  below  is  essentially  the  same  as  in  [8]  ex- 
cept that  in  [8]  d is  assumed  to  be  two  while  here  it  can  be  any  positive  integer. 

The  algorithm  is  applied  to  a parametric  form  of  an  n by  dn  LCP: 

d 


(3.13a) 

(3.13b) 

(3.13c) 


w = q + Xp  + £ M X 
i = l 


w > 0,  0 < xi  < r*,  4i  = l,...,d-l,  x > o 


wV  = 0,  (r*  - x*)T  xi  + 1 = 0, 


i = 1, . . . ,d-l  , 


where  X c (0,  X),  X > 0 • 

0 d 

In  a general  iteration  step,  there  exist  index  sets  y ,...,Y  such  that 


basic 


w . = 0:  nonbasic, 
Y3 


and  for  i = 1 , . . . ,d 


x1  =0:  nonbasic, 
Y3 

x1  : basic 


j = 1, . . . ,d 


j = 0, . . . ,i-l 


x1  = r1  . : nonbasic,  j = i+l,...,d 


We  shall  denote  by  q,p,M,  the  current  values  of  the  data.  Also  we  denote  by 
the  n-vector  defined  by 

d-1 


= l l . r . M1  . 
i=l  jcy1  3 ° 


Algorithm  II-  Condensed  Graves'  algorithm. 

Step  1:  If  P Q >0,  P d > 0 and  p . = 0,  j = 1,.. 

y Y Y 

basis  gives  the  solution  at  X = X . 

Step  2:  Determine  the  next  critical  index  s by 
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. , d- 1 , then  stop:  the  currc 


1 


U*tk 


f 


1 


P3  < 0,  j - 1, ...  ,11  ; 


s c arg  min  / 


r . - (q. +y . ) ii  1 

_2 X p > 0,  j £ <J  v1 

3 ] i=l 


The  next  critical  value.  X',  is  given  bj 


it  p < 0 

s 

if  p > 0 

s 


Step  3 : If  X*  > X,  then  stop:  the  current  basis  cjives  the  solution  at  X - X 

Step  4 : Let  s be  the  critical  index  determined  in  Step  2.  The  next  pivot  element  is 

* 

i * 

Mg  t where  i is  given  by  the  following: 

i * 

(i)  p < 0,  s c y , i c 10,1},  then  i = 1 

i * 

(ii)  p < 0,  c c y , r » (2,.  , d},  then  i - i-1 

I * 

(iii)  ps  > 0,  s e Y , it  {!,..., d-1),  then  i = iPl  . 

Determine  the  new  "eta-vector”  by  computing  the  current  value  of  M1  ("forward 
transformation").  Update  q,  p and  y by  standard  pivot  formulas.  Heturn  to  Step  1. 


3.4.  Comparison  of  Algorithms.  Suppose  that  we  are  to  solvi  th'-  problem  of  determining 
the  behavior  of  a structure  with  n critical  points  of  which  m n ■ "redundant”  and 
hj  (Xj)  having  d >1  linear  pieces,  j l,...,n  . Our  objective  is  to  compare  the 
computer  storage  space  and  the  number  of  operations  required  to  solve  the  problem  by 
(i)  hemkc ' s algorithm  (on  (1.0)),  (ii)  the  compact  inverse  algorrttim  (on  (3.8))  or 
liii)  tlie  condensed  Graves'  algorithm  (on  (3.12)). 

The  results  of  the  compar  ison  are  summarized  in  Tabic  -1.1.  To  make  the  comparison, 
it  is  as  imccl  that  specifications  of  the  first  two  algorithms  arc  as  stated  in  Section  2, 
and  those  of  the  thrrd  one  are  the  same  as  Lcmkes'. 

It  is  clear  that  for  d > 1 the  condensed  Graves'  algorithm  is  always  superior  to 

terokes'  with  respect  to  both  criteria,  provided  that  T'  and  T are  close”  Thus, 

0 Refer  to  footnote  1. 
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■ I.I.Ch. 


Lemke 


compact  inverse 


condensed  Graves 


storage 


„2  2 

d n + T dn  + dn 


dmn  + 3dn  + 


m (m+I ) 


n + Tn  + 2dn 


operation 

counts 


T' (T'-l) 

dn  -L-~ — +3T'dn 


T dn (m+3) +2  (mi  1) 


T(T-l) 

n — — + 4Tn 


Table  3.1.  Storage  and  operation  counts  for  three  algoritlims 

Lemke's  algorithm  will  be  omitted  from  further  considerations. 

The  relative  advantage  of  the  two  special-purpose  algorithms  depend  on  the  values  of 
n,m,d  and  T.  We  shall  discuss,  below,  likely  values  of  n,m,d  and  T in  practical  problems, 
but  in  general  the  following  may  be  asserted:  The  condensed  Graves'  algorithm  is  more 
advantageous  when  m/n  is  large  and  T is  small,  and  the  compact  inverse  algorithm  is 
more  advantageous  when  m/n  is  small  and  T is  large. 

In  this  particular  application,  the  matrices  involved  (in  (3.8)  and  (3.12))  tend  to 
bo  highly  dense. 

The  value  of  n could  be  quite  large  if  for  instance  the  underlying  structure  is 
continuum  In  practical  reinforced  concrete  frame  problems,  small  to  medium  size  n is 
expected. 

The  value  of  m varies  from  one  problem  to  another,  but  the  ratio  m/n  is  expected 
to  be  in  the  neighborhood  of  one  half  in  most  cases. 

The  value  of  d depends  on  how  accurately  the  piecewise  linear  approximation  of  h 
must  be  done.  In  reinforced  concrete  frame  problems,  for  example,  d = 2 or  3 seems  ade- 
quate . 

Finally,  in  this  application  the  number  of  pivots,  T,  is  expected  to  be  of  a moder- 
ate magnitude,  even  when  the  parametric  problem  (3.7)  must  be  solved,  the  range  of  X is 
relatively  small.  The  actual  value  of  T could  be  highly  variable  but  in  general  it  is 
expected  to  be  large  when  dn  is  large  (note  that  in  both  (3.8)  and  (3.12)  there  are  dn 
variables) 
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3.5.  Computational  Result:  Based  on  the  comparison  in  terms  of  the  operation  counts, 

it  is  expected  that  the  condensed  Graves'  algoritlim  solves  the  problem  more  quickly  than 
the  compact  inverse  algorithm,  when  m/n  is  large  and  the  relative  advantage  is  reversed 
when  m/n  l ■•comes  small  (such  a value  of  m/n  will  be  called  the  reversal  value)  . 

In  order  to  (i)  verify  this  tendency  and  (ii)  determine  the  reversal  value  of  m/n  , 
we  have  i i>r  1 < u mod  a compr.  ational  experiment  comparing  l lie  computational  times  required  by 
the  two  algor i thins  to  solve  some  problems  of  various  sizes.  The  results  are  reported 
below. 

We  fi.it  c.  ..  .tructvd  and  solved  the  n by  dn  LCP  (3.12)  by  the  condensed  Graves' 
algorithm  lot  t ..  : et:  at  data  with  various  values  of  d,  n and  m (see  below).  The  data 
were  randomly  . ■ i.ited  in  the  following  form  with  the  specified  properties: 

q > 0,  p = -NsE  < 0 , 

i i T 

M = U + GG  with  n by  m G and 
1 i i i. 

II  diaqt  j , . . . ,u  1,  l l,...,d,  where 
1 n 

1 d 

Uj  > . . . > u_.  >0,  j = 1, . . . ,n 

rX  > 0,  i = 1, . . . ,d-l  . 

We  then  co;.  cited  the  data  to  obtain  the  corresponding  LCP  (3.8),  and  solved  it  by 
the  compact  inverse  algc.i  ithm. 

A total  of  147  problems  were  solved  by  each  of  the  two  algorithms  with  the  follow- 
ing values  of  d,n  and  m;  problems  with  the  same  values  of  d and  n are  considered  to 
be  in  one  grout 


Group 

1: 

d = 

2, 

n 

« 40, 

m = 

2,4,28,30 

and 

m 

= 

8 

+ 

2t, 

t 

0,. 

. . ,6 

Group 

2: 

d = 

2, 

n 

= 50, 

m = 

2,4,38,40 

and 

m 

= 

8 

+ 

2t, 

t 

= 0, 

. ..,10 

Group 

3: 

d = 

2, 

n 

= 60, 

m - 

2,4,48,50 

and 

m 

= 

8 

2t , 

t 

= o. 

• ••,13 

Group 

4: 

d - 

2, 

n 

70, 

n = 

2,4,48,50 

and 

m 

= 

10  t 

2t , 

t 

= 0, 

. ..,15 

Group 

5: 

d = 

3, 

n 

10, 

m e 

2,4,28,30 

and 

m 

= 

8 

2t , 

t 

■=  0, 

Group 

6: 

d 

3, 

n 

- 50, 

m 3 

2,4,38,40 

and 

m 

9 

+ 

3t. 

t 

= 0, 

...,7 

Group 

7: 

d - 

3 , 

n 

60, 

m * 

2,4,48,50 

and 

m 

= 

9 

+ 

3t , 

t 

= 0, 

• • • ,7 

Croup 

8: 

d 

3, 

n 

= 70, 

m - 

2,4,48,50 

and 

m 

4= 

9 

4 

3t, 

t 

= 0, 

...,9 

We  wrote  a Fortran  program  for  the  condensed  Graves'  algorithm  which  was  based  on 
the  specifications  stated  in  Section  4.3.  The  program  for  the  com,  act  inverse  algorithm 
used  was  essentially  the  same  as  that  mentioned  in  Section  3,  except  for  a certain  modi- 
fication which  takes  into  account  the  special  forms  of  the  matrices  P and  N . 

Results  01  the  experiment  are,  in  general,  closely  in  accordance  with  the  analysis 
based  on  the  operation  counts.  Some  of  them  are  shown  and  explained  in  the  following. 

Figures,  plotting  the  computation  times  (exclusive  of  input/output ) per  iteration 
in  both  algorithms  are  diawn  for  problems  in  each  of  the  eight  groups.  Figures  3.1- 
3.3  are  those  for  Groups  1,  4 and  7;  others  are  omitted  since  they  appear  quite  similar  to 
those  shown  hire.  The  above-mentioned  tendency  of  the  relative  advantage  of  the  algorithms 
is  clearly  demonstrated  in  all  the  figures. 

These  figures  also  confirms  a result  of  the  analysis  in  Table  3.1  that  the  time  per 
iteration  in  the  compact  inverse  algorithm  increases  quadratical ly  in  m . Figure  3.4 
shows  that  the  time  per  iteration  in  the  condensed  Graves'  algorithm  increases  linearly 
in  T,  confirming  the  corresi'onding  analytical  result  in  Table  .'.1 

In  Figures  3. 1-3. 3 the  time  per  iteration  in  the  condensed  Graves'  algorithm  appears 
to  be  decreasing  in  m and  thus  in  conflict  with  the  analytical  result  (which  is  inde- 
pendent of  m ).  This  was,  however,  a natural  consequence  of  the  fact  that  T,  the 
number  of  iterations  needed  to  solve  the  problem,  tended  to  be  smaller  as  m became  larger 
recall  that  the  time  per  iteration  increases  linearly  ir,  T . 

Approximate  values  of  reversal  were  determined  in  the  following  fashion.  In  Figures 
3. 1-3. 3 (and  in  similar  ones  not  shown  here)  we  drew  horizontal  lines  to  indicate  the  time 
per  iteration  in  the  condensed  Graves'  algorithm  corresponding  to  a problem  which  is  solved 
in  T iterations,  where  T is  the  average  of  the  numbers  of  iterations  necessary  to 
solve  the  problems  in  Group  1,  4 and  7,  respectively.  We  need  to  do  so  since  the  time  per 
iteration  in  the  condensed  Graves'  algoritlim  depends  linearly  in  T (as  shown  in  Figure 
3.4)  and  the  number  of  iterations  needed  to  solve  each  of  the  problems  in  each  group 
varied. 

7 This  tendency  may  be  attributed  to  the  manner  in  which  the  data  were  constructed;  we 
do  not  elaborate  on  this  point  here,  however. 


In  each  of  Figures  3. 1-3.  3 let 


m be  the  value  of  m corre;  ponding  to  the  "cross- 


over" point,  i.e.  the  point  at  which  the  horizontal  line  meet-  the  quadratic  curve 

corresponding  to  the  time  for  compact  inverse  algoritlim.  For  the  reason;,  mentioned  above, 

* 

it  is  appropriate  to  consider  m /n  as  an  approximate  value  of  reversal  Table  3.2  lists 
approximate  values  of  reversal  obtained  in  this  fashion  along  v . th  other  ir.b  rit.it  ion  for 
problems  in  each  group. 


Group 

i 

2 

3 

4 

5 

0 

7 

8 

reversal 

* 

, m 

value  — 
n 

.3 

. 3 

3 

. H 

.4 

. S 

.4 

•4 

average  # 
of  iterations 

38 

SI 

70 

82 

71 

100 

83 

83 

d 

2 

2 

2 

2 

3 

3 

3 

3 

n 

40 

so 

60 

70 

40 

. 

so 

GO 

70 

Table  3.2.  Approximate  reversal  values. 

Wc  parenthetically  mention  that  the  reversal  value  tends  to  be  large  when  the  average 
n umbct  of  iterations  is  large;  this  is  consistent  to  the  analytical  result  in  Table  3.1. 

For  further  computational  results,  we  refer  to  a report  [101  by  one  of  the  authors 
establishing  the  predicted  superiority  of  the  condensed  Graves'  algorithm  over  benthos' 
with  respect  to  the  computation  lime. 

3.(>  Cone -In-. ion;  . In  the  structural  engineering  application,  we  can  draw  the  following 

conclusions  on  the  suitability  of  the  solution  methods. 

(i)  The  condensed  Graves'  algorithm  is  always  better  than  r.etnkt  's  algorithm,  both 
in  terms  of  the  computer  storage  and  computation  time  required  to  solve  the  problem. 

(ii)  In  most  problems  (where  m is  about  n/2  ),  the  condensed  Graves'  algoritlim 
solve  the  problem  more  quickly  requiring  more  storage  space  than  the  compact  inverse  al- 
gorithm. 

(iii)  The  compact  inverse  algorithm  is  better  than  the  condensed  Graves'  algorithm  in 
terms  of  both  storage  and  computer  time  r£  m/n  is  small  and  dn  j large. 
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4.  GRADUATION  APPLICATION 


4.1.  Problem  Description.  In  this  section,  we  discuss  the  third  and  final  application  of 
the  class  of  LCP's  and  PLCP's  with  matrices  of  the  form  (1.1).  This  application  arisen  in 
the  area  of  graduation  or  curve  fitting  and  smoothing.  Before  describing  how  the  comple- 
mentarity problems  emerge  we  explain  the  meaning  and  process  of  graduation  in  the  context 
of  actuarial  sciences.  The  following  discussion  is  based  largely  on  the  excellent  treatise 
of  the  subject,  a monograph  by  M.  D.  Miller  (18). 

The  actuary  is  concerned  with  the  contingencies  of  death,  disability,  retirement, 
sickness,  withdrawal,  marriage,  etc.  He  must  know  the  probabilities  of  such  events  in 
order  to  predict  their  future  occurrence  and  in  order  to  be  able  to  calculate  premiums,  re- 
serves, annuities  and  so  forth,  for  insurance  and  other  financial  institutions.  Tables 
setting  forth  such  probabilities  must  be  constructed  and  for  that  purpose,  observations 
are  made  of  the  happening  of  such  events.  Graduation  is  one  of  the  steps  in  the  construc- 
tion of  these  tables. 

Typically,  a series  of  observed  probabilities  is  found  to  contain  irregularities.  It 
is  with  irregular  series  of  observed  values  of  continuous  varying  quantities  such  as  mor- 
tality rates,  that  the  problem  of  graduation  deals.  Formally,  graduation  may  be  defined  as 
the  process  of  securing  from  an  irregular  series  of  observed  values  of  a continuous  vari- 

1 

able  a smooth,  regular  series  of  values  consistent  in  a general  way  with  the  observed 
series  of  values.  This  smooth  series  of  graduated  values  is  then  taken  as  a representation 
of  the  underlying  law  which  gave  rise  to  the  series  of  observed  values. 

There  are  several  theoretical  reasons  justifying  the  postulate  of  a smooth  regular 
underlying  series  of  true  values  in  the  process  of  graduation.  One  of  these  reasons  can 
be  dealt  with  using  the  mathematical  theory  of  probability,  see  (18]. 

Among  those  methods  by  which  the  process  of  graduation  may  be  accomplished  is  the 
difference-equation  method.  This  method  is  founded  on  the  formulation  of  an  analytic  ex- 
pression measuring  the  combination  of  two  essential  characteristics  of  graduation: 

(i)  smoothness  and  (ii)  fitness,  or  consistency  with  the  observed  data.  On  the  one  hand, 
the  graduated  series  should  be  smooth  compared  with  the  unqraduated  series;  on  the  other 
hand  it  should  be  consistent  with  the  indications  of  the  ungraduated  series.  Basically, 


I 
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these  two  qualities  ol  smoothness  and  fitness  are  inconsistent  in  the  sense  that  smooth- 


ness may  not  he  improved  without  some  sacrifice  of  fitness,  and  vice  versa.  This  situation 
is  similar  to  Markowitz'  mean-variance  model  in  portfolio  analysis. 

Numerically,  smoothness  and  fitness  are  measured,  respectively,  by 


(4.1) 


F = £ w (u  - u")“ 


x=l 


and 


(4.2) 


n-m 

s = y (Am  u )■ 
x=l  X 


whei  o 


u£  = ungraduated  values  for  integral  values  x = 1 to  x = n 
u = con  t-spond  ing  graduated  value 
w = weight  assigned  to  the  value  u” 


and 


A u - m-th  forward  difference  of  the  graduated  series 
x 

\ 

m 

= l (-l)m'r  (m)u  ^ 
r 0 1 X+r 


m integer  indicating  the  desired  degree  of  smoothness;  typical  values  are 
?,  3,  4 arid  5 . 


If  the  ungraduated  values  u^  represent  the  mortality  rates,  then  the  weights  w^  are 
normally  assigned  the  corresponding  numbers  exposed  to  the  risk  of  death. 

For  a specified  value  of  in,  the  difference-equation  graduation  formula,  also  re- 
ferred to  as  Win t ♦ aket -Henderson  formula,  is  derived  from  the.  minimization  of  the  expres- 
sion F ♦ kS  where  k is  the  constant  indicating  the  relative  emphasis  we  are  placing  on 
fitness  and  .mo  thner.'  Before  the  advent  of  electronic  computers,  it  was  generally  felt 
that  this  grab  tion  method,  while  producing  excellent  results,  involved  a working  process 
so  long  and  difficult  \.  to  be  impra<  tical  except  for  making  graduations  of  considerable 
importance.  Nevertlu  less,  as  it  has  been  pointed  out  in  (7),  this  method  is  really  very 
simple,  especially  when  v>.  tor-matrix  notation  is  employed.  Indeed,  letting  u"  and  u 


be  the  vectors  of  ungraduated  and  graduated  values  respectively,  1 be  the  diagonal  matrix 

m 

of  weights,  and  G be  the  coefficient  matrix  of  the  m-th  order  difference  operator  A , 

i.e. 


(4.3)  G - 


we  have 

(4.4) 


<-l)m 

0 

m. 

. m 

(-1) 

<x)  (-D 

m-  2 

m. 

. . m- 1 

(-1) 

v (~1> 

0 

(-i)m  0 

<-nm_1  <■)  (-Dm 


T T T T 

F ( ki  u (1  ♦ kGG  ) u - 2u  i.  u"  + u"  u" 


Differentiation  with  respect  to  u gives 


(Y  + kGG  )u  ■=  l u" 


from  which  the  graduated  series  (u^l  can  be  obtained  readily  Observe  that  G is  of 
order  n by  (ri-m)  and  the  nonzero  elements  of  each  of  its  column  are  defined  by  the  (m  + 1) 
binomial  coefficients  (-1)  , (-1)  ' (^1,  •••*  ^ • Moreover,  each  row  of  G contains  at 
most  (in  + 1)  non-zero  entries. 

Summing  up  the  above  discussion,  we  see  that  the  difference-equation  approach  for 
graduation  involves  nothing  but  the  solution  of  a system  of  linear  equations  (4.S)  derived 
from  an  unconstrained  minimization  piohlem. 


Despite  the  fact  that  this  particular  appi.-.u  h for  «u  * i •- 1 • . it  j i..t  i • • i.  widely  ac- 
cepted and  used  by  researchers  and  pra-  ' iticnci  in  actuari  il  < • , we  ft  . i that  it  is 

imperfect  and  could  be  improved.  In  particular,  w«.-  proj  > •.  « t!<  t.  a t » .trained  dif- 

ference-equation approach:  Let  C denote  t he-  set  of  a p.  i »i  1 \i,u  it  i-  that  v-  hav- 
about  the  series  of  ungraduated  values  , i.e.  C is  a give.,  i t t Rn  antaini  r*g 
the  vector  u"  . The  new  approach  requires  that  the  graduated  lx  *.  u be  th*  sc  lutiun 

to  the  following  constrained  minimization  problem 

(4.6)  minimize  F t kS  subject  to  u e C 


Tin.  reason  for  consideri ng  uch  a constrained  problem  is  cb  a.  . . j . i ) •-.••ate  ...  feel 

that  the  graduated  series  should,  in  addition  to  exhibiting  con  tency  and  tl  • a me 

tendency  as  the  ungraduated  series,  possess  tin  same  characteristics  as  the  latter  series. 
For  instance,  if  vs  are  graduating  mortality  rates,  we  certainly  hope  that  the  graduated 
valuer,  v.-  ill  roji  os<  rit  probabilities.  By  taking  C to  Lh  the  unit  cub*  in  (4.6),  we  are 
guarani  ed  th  it  t la  graduated  valu*  wi  1 1 indeed  )•«  • nonnegati  ml  . * gre  itei  than  or..-. 
It  should  be  j>oint"d  out  that  requiring  th  graduated  values  to  be  pr* »haLi  1 i tics  by  no 
mean  i mi  lies  that  the  graduated  series  is  a probability  den  ity  fun-ton.  Therefore,  no 
restriction  on  the  sum  of  graduated  values  is  necessary. 

In  the  ing  irtant  case  where  C is  a rectangle,  we  can,  by  an  c . , < -Lange  of  vai l - 
abU  .,  cast  problem  (4.6)  in  t.h<  form  below: 


(4.7) 


...IT  T T 

minimize  u (>.  + kGG  )u  + q u subject  to 


0 < u < a 


To  the  authors’  knowledge,  problem  (4.7)  has  not  explicitly  apj  ned  in  the  literature 
Nevertheless,  several  of  its  special  cases  have  been  studied  under  various  contexts.  Of 
course,  without  the  constraints,  problem  (4.7)  reduces  to  th  classical  difference  equation 
approach.  Incidentally,  the  unconstrained  problem  has  also  1 * en  studied  in  the  context  of 
spline  function'  , see  U3J,  (20),  (21),  (23).  Mangasarian  and  . •humak*  i 114)  have  pro- 
vided necessary  and  sufficient  optimality  conditions  for  the  ca a ) 0 in  terms  of 

difference  operators.  The  continuous  version  of  this  special  can*  hu*  appeared  in  (IS) 
under  the  context  of  control  theory  and  in  (3)  under  the  context  f mi  engineering  appli- 
cation . 
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The  Karush-Kuhn-Tucker  optimality  conditions  for  (4.7)  are  given  by: 

(4.8a)  w = q + (5)  4 kGGT) u + v >_  0 , u _>  0 

(4.8b)  z = a -u  2.0,  v j>  0 

(4.8c)  v u = z^v  - 0 

If  the  upper  bounds  a^  are  all  infinity,  then  problem  (4.8)  reduces  to  a linear  comple- 
mentarity problem  (q,M)  where  M is  of  the  form  (1.1).  Therefore  the  compact  inverse 
algoritlim  is  directly  applicable  in  this  case.  The  general  case  where  not  all  the  a 's 
are  infinity  can  be  treated  by  an  extension  of  the  algoritlim.  See  [19]. 

T 

In  what  follows,  we  point  out  some  characteristics  of  the  matrix  £ + kGG  appearing 
in  this  graduation  application  which  can  be  used  advantageously  in  the  application  of 
Algorithm  I (and  its  extension).  Since  m is  typically  very  small,  the  n by  (n-m)  matrix  is 
almost  square.  According  to  the  comparison  table  in  Section  2.3,  this  latter  fact  would 
seem  to  imply  that  the  application  of  the  compact  inverse  algorithm  would  bo  undesirable. 
Nevertheless,  (4.3)  shows  that  G is  very  sparse  and  specially  structured.  Indeed,  it  ir. 

easy  to  see  that  for  every  B in  {l,...,n},  the  matrix  A(B)  defined  in  (2.9)  is  a 

* 

band  matrix  with  band  width  2m  4 1 Therefore  the  Cholesky  factor  L has  at  most  m 

nonzero  off-diagonal  elements  in  each  row.  The  sparsity  of  L reduces  the  amount  of  op- 

2 

erations  in  Subroutines  I and  II  from  2 (n-m)  + 4 (n-m)  to  n(4m+f>)  - 7m(m+l)  . Moreover, 

the  fact  that  the  nonzero  elements  in  each  row  of  G are  coefficients  of  the  m-th  ordci 
m P1  k 

difference  operator  A makes  it  clear  that  the  vector  £,  ' defined  by  (2.22)  can  be 
calculated  without  any  multiplication  Finally,  because  of  its  -special  structure,  no 
storage  space  is  required  for  the  matrix  G in  the  implementation  of  the  algorithm.  Its 

elements  can  be  generated  very  easily  if  needed. 

T 

The  matrix  M = £ + kGG  is  sparse  if  G is  (although  M is  less  sparse  than 
G)  , and  this  could  certainly  be  taken  into  account  to  increase  the  efficiency  of  general- 
purpose  algorithms  such  as  Lomke's.  However,  the  fact  that  the  compact  inverse  algorithm 
could  be  implemented  so  that  it  takes  full  advantage  of  the  special  structures  of  G in- 
dicates  the  possibility  that  the  special-purpose  algoritlim  may  turn  out  to  be  more  effi- 
cient than  general-purpose  algorithms.  In  order  to  draw  empirical  conclusions,  we  have  per- 


formed several  computat ional  experiments.  The  results  are  reported  in  the  next  subsection. 
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4.2.  Computat  ion  Results.  In  order  to  compare  the  efficiencies  of  the  compact  inverse 
and  Lemke's  algorithms  in  solving  (P)  LCP’s  with  matrix  of  the  form  (1.1),  where  G is 
given  explicitly  by  (4.3),  experiments  are  designed  to  solve  the  following  variant  of 
problem  (4.7) : 

IT  T T 

(4.9)  minimize  j u (1.  + KGG  )u  + q u subject  to  u _>  0 . 

We  do  not  consider  upper  bounds  in  the  experiments  because  the  code  NULEMKE  that  we  use 
for  Lemke's  algorithm  handles  upper  bounds  as  separate  constraints,  whereas  the  compact 
inverse  algorithm  treats  them  implicitly.  It  is  certainly  disadvantageous  to  Lemke's  al- 
gorithm if  the  upper  bounds  are  included.  It  should  be  mentioned  that  NULEMKE  does  take 
advanLage  of  sparsity  of  the  data. 

A set  of  three  experiments  corresponding  to  m = 2,3  and  4 are  performed.  For  each 

8 

value  of  m,  five  values  50,  60,  70,  80  and  90  of  n are  tested.  In  the  case  m = 3 , 
it  is  found  that,  if  not  for  the  reinversion  in  Lemke's  algorithm,  which  occurs  when 
n = 90,  the  two  algorithms  produce  results  very  close  to  each  other.  The  results  for 
m = 2 and  4 are  summarized  in  Tables  4.1  and  4.2. 


n 

density  of 

M (%) 

(1  of 

iterations 

2 -3 

time/iteration  (10  sec.) 

Lemke  compact  inverse 

50 

9.8 

48 

11.6 

14.0 

60 

8.2 

58 

15.) 

16.7 

70 

7.0 

67 

16.7 

19.3 

80 

6.2 

76 

18.9 

21.9 

90 

5.5 

88 

29. 33 

24.7 

Table  4.1  m = 2 


8 These  values  of  m and  n are  typical  in  actuarial  graduation. 


n 

density  of 

M (%) 

# of 

iterations 

time/iterati 

Lemke 

2 

on  (10  sec.) 

compact  inverse 

50 

17.2 

58 

18.1 

18.6 

60 

14.4 

70 

26.0 

23.0 

70 

12.5 

B4 

27.8 

25.4 

80 

10.9 

109 

36. 73 

29.6 

90 

9.7 

114 

38. 23 

33.3 

Table  4.2  m = 4 


Notes:  l.We  chose  the  covering  vector  in  Lemke's  algorithm  so  that  the  numbers  of  itera- 
tions in  both  algorithms  are  the  same. 

2.  The  time  is  exclusive  of  input/output. 

3.  The  sudden  increase  is  due  to  the  occurrence  of  reinversion  of  the  basis. 

‘5.3.  Conclusions.  Based  on  the  above  computational  results,  we  may  draw  the  following 
three  conclusions: 

(i)  The  compact  inverse  algorithm  requires  much  less  storage  space  than  Lemke's 
algorithm. 

(ii)  Although  the  ratio  (n-m)/n  is  close  to  one  in  these  graduation  problems,  the 
compact  inverse  algorithm  is  performjng  compatibly  with  Lemke's  algorithm,  in  terms  of 
computation  time. 

(iii)  The  results  indicate  that  as  m gets  larger,  or  equivalently,  as  G becomes 
"thinner",  the  efficiency  of  the  compact  inverse  algorithm  becomes  higher,  a fact  consistent 
with  the  analysis  of  the  algorithm  made  in  Section  2.3. 
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